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Abstract 

In this paper we formulate Bulgac-Kusnezov constant temperature dynamics in phase space 
by means of non-Hamiltonian brackets. Two generalized versions of the dynamics are similarly 
defined: one where the Bulgac-Kusnezov demons are globally controlled by means of a single 
additional Nose variable, and another where each demon is coupled to an independent Nose-Hoover 
thermostat. Numerically stable and efficient measure-preserving time-reversible algorithms are 
derived in a systematic way for each case. The chaotic properties of the different phase space flows 
are numerically illustrated through the paradigmatic example of the one-dimensional harmonic 
oscillator. It is found that, while the simple Bulgac-Kusnezov thermostat is apparently not ergodic, 
both of the Nose-Hoover controlled dynamics sample the canonical distribution correctly. 
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I. INTRODUCTION 



In condensed matter studies, there are many situations in which molecular dynamics sim- 
ulation at constant-temperature [l-3] is needed. For example, this occurs when magnetic 
systems are modelled in terms of classical spins {4- 7]. Deterministic methods 0-0, based 



on non-Hamiltonian dynamics |llHl9|. can sample the canonical distribution provided that 
the motion in the phase space of the relevant degrees of freedom is ergodic jl, 3]. H< 



classical spin systems are usually formulated in terms of non-canonical variables 



owever, 



20|,|2l|, 



without a kinetic energy expressed through momenta in phase space, so that Nose dynamics 



cannot be applied directly. To tackle this problem, Bu 



22 



gac and Kusnezov (BK) introduced 
24| which can be applied to spins. A 



a deterministic constant-temperature dynamics 
number of numerical approaches to integration of spin dynamics can be found in the lit- 
erature 



25 



28| . However, BK dynamics, as any other deterministic canonical phase space 



flow, is able to correctly sample the canonical distribution only if the motion in phase space 
is ergodic on the timescale of the simulation. In general, this condition is very difficult to 
check for statistical systems with many degrees of freedom, while it is known that, despite 
its simplicity, the one-dimensional harmonic oscillator provides a difficult and important 
challenge for deterministic thermostatting methods {9! [29^ 31 ] . 

In this paper, we accomplish two goals. First, by reformulating BK dynamics through 
non-Hamiltonian brackets Q, Q in phase space, we introduce two generalized versions 
of the BK time evolution which are able to sample the canonical distribution for a stiff 
harmonic system. Second, using a recently introduced approach based on the geometry 



of non-Hamiltonian phase space 19], we are able to derive stable and efficient measure- 
preserving and time-reversible algorithms in a systematic way for all the phase space flows 
treated here. 

The BK phase space flow introduces temperature control by means of fictitious coor- 
dinates (and their associated momenta in an extended phase space) traditionally called 
'demons'. Our generalizations of the BK dynamics are obtained by controlling the BK 
demons themselves by means of additional Nose- type variables jsj|. In one case, the BK 
demons are controlled globally by means of a single additional Nose-Hoover thermostat 8|,19]. 
In the following this will be referred to as (BKNH) Bulgac-Kusnezov-Nose-Hoover dynam- 
ics. In the second case, each demon is coupled to an independent Nose- Hoover thermostat. 
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This will be called the Bulgac-Kusnezov- Nose- Hoover chain (BKNHC), and corresponds to 
'massive' NH thermostatting of the demon variables j^ . 

The ability to derive numerically stable measure-preserving time-reversible algorithms 



19[ for Nose controlled BK dynamics is very encouranging for future applications to ther- 



mostatted spin systems. 

This paper is organized as follows. In Sec. HI] we briefly sketch the unified formalism for 
non-Hamiltonian phase space flows and measure-preserving integration. The BK dynamics 
is formulated in phase space and a measure-preserving integration algorithm is derived in 
Sec. IIIIL The BKNH and BKNH-chain thermostats are treated in Sees. [IV] and Irrespectively. 
Numerical results for the one-dimensional harmonic oscillator using these thermostats are 
presented and discussed in Sec. IVT] Section IVHI reports our conclusions. 

In addition we include several Appendices. A useful operator formula is derived in Ap- 
pendix [A], while invariant measures for the BK, BKNH, and BKNHC phase space flows are 
derived in Appendices EJ EH and U2 respectively. 
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II. NON-HAMILTONIAN BRACKETS AND MEASURE-PRESERVING ALGO- 
RITHMS 



Consider an arbitrary system admitting a time-independent (extended) Hamiltonian ex- 



pressed in terms of the phase space coordinates Xi 



1 2N. In this case, the Hamil- 



21 



tonian can be interpreted as the conserved energy of the system. 

Upon introducing an antisymmetric tensor field (generalized Poisson tensor 
phase space, B(x) = —B T (x), one can define non-Hamiltonian brackets 14j-|l6| as 

r za sr^ oa db 



33|) in 



(1) 



where a = a(x) and b = b(x) are two arbitrary pha se space functions. The bracket defined 
in Eq. ([I]) is classified as non-Hamiltonian 14j-ll6| since, in general, it does not obey the 



Jacobi relation, i.e., in general the Jacobiator J ^ 0, where 21] 



J = {a, {b, c}} + {b, {c, a}} + {c, {a, b}} 



(2) 



with c = c(x) arbitrary phase space function (in addition to the functions a and b, previo usly 
introduced). If J ^ 0, the tensor Bij is said to define an 'almost-Poisson' structure [34]. 
(Such systems have also been called 'pseudo-Hamiltonian' |33|.) 

An energy- conserving and in general non-Hamiltonian phase space flow is then defined 
by the vector field 

(3) 



3=1 



where conservation of H(x) follows directly from the antisymmetry of B^. 

It has previously been shown how equilibrium statistical mechanics can be comprehen- 
sively formulated within this framework 16fl . It is also possible to recast the above formalism 

nn 

and the corresponding statistical mechanics in the language of differential forms [121, [18] . If 
the matrix B is invert ible (this is true for all the cases considered here), with inverse fi^, 



we can define the 2-form 



35] 



Q = \Vt i:j dx L A dx j . 



(4) 



The dynamics of Eq. ([3]) is then Hamiltonian if and only if the form fll]) is closed, i.e., has 



zero exterior derivative, dQ = 



351 ] . This condition is independent of the particular system 



of coordinates used to describe the dynamics. 



The structure of Eq. can be taken as the starting point for derivation of efficient 
time-reversible integration algorithms that also preserve the appropriate measure in phase 
space [3]. Measure-preserving algorithms can be derived upon introducing a splitting of 
the Hamiltonian 

n s 

H = J^H a (5) 

0=1 

which in turn induces a splitting of the Liouville operator associated with the non- 
Hamiltonian bracket in Eq. (JTJ) , 

2N dH 

L a Xi = {xj, H a } = ^ ] Bjj ^ ■ (6) 
j=i 1 

When the phase space flow has a non-zero compressibility 

4-^ dxi dxn 

1,3=1 J 

;he statistical mechanics must be formulated in terms of a modified phase space measure 
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-lis I 

uJ = e ~ w{x) u (8) 

where 

u = dx 1 A dx 2 A ... A dx 2N (9) 

is the standard phase space volume element (volume form [35() and the statistical weight 
w(x) is defined by 

dw , . , . 

- = «(,) . (10) 

It has been shown that, provided the condition 

A =0, z = l,...2iV (11) 

is satisfied, then 

L Q cJ = for every a , (12) 

I — I 

so that the volume element ZJ is invariant under each of the L a [19]. The condition ( II ip is 
satisfied for all the cases considered below, so that, exploiting the decomposition in Eq. (JOJ), 
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algorithms derived by means of a symmetric Trotter factorization of the Liouville propagator: 



n 3 — 1 



a=l 
n 3 — 1 

x ] [ ex P 

0=1 



exp[rL] = JJ exp expexp[rL„J 



(13) 



are not only time-reversible but also measure-preserving. 
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III. PHASE SPACE FORMULATION OF THE BK THERMOSTAT 



A phase space formulation of the BK thermostat can be achieved upon introducing the 
Hamiltonian 

KM) , K 2 {p^) 



H 



BK 



f - + V(q) + 
2m 



+ 



k B T(( + 



(14a) 
(14b) 



where (q,p) are the physical degrees of freedom (coordinates and momenta), with mass 
m, to be simulated at constant temperature T, while ( and £ are the BK 'demons', with 



corresponding inertial parameters and m^, and associated momenta {p^,p{) |22|-|24j. K\ 
and K 2 provide the kinetic energy of demon variables, and for the moment are left arbitrary. 

Upon defining the phase space point as x = (q, £, C,p,Pc,p$) = x 2 , X3, %4, %e) , one 
can introduce an antisymmetric BK tensor field as 



B 



BK 








-1 



G 2 








dGi 
dp 

- 









dG 2 
dq 



1 





d 






dGi 
dp 








-G 2 


dG 2 
dq 







(15) 



where G\ and G 2 are functions of system variables (p, q) only. 

Substituting B BK and H BK into Eq. fl3]), we obtain the energy-conserving equations 

dH G 2 (q,p)dK 2 
dp 



Q 

C 

£ 

V 



dp£ 
1 dGidKx 
rriQ dp dp^ 
1 dG 2 dK 2 
dq dp^ 
dH Gi(q,p) dK x 
dq dpc_ 
dH dGi 

G l"a k ^ T ~^— 

dp dp 



G 



dH 



dG 2 



dq " dq 

The associated invariant measure for the BK flow is discussed in Appendix IB1 



(16a) 
(16b) 
(16c) 
(16d) 
(16e) 
(16f) 
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A. Algorithm for BK Dynamics 



In order to derive a measure preserving algorithms, the first step, following Eq. (jSJ), is to 
introduce a splitting of H BK : 



Hf K = V(q) (17a) 

H* K = (17b) 

Hf K = k B T( (17c) 

Hf K = k B Ti (17d) 

^BK = KM (1?e) 

Hf K = . (17f) 

A measure-preserving splitting of the Liouville operator then follows from Eq. ([6]): 

1 dq dp dq dp^ 

lf = 4 + C^ (18b) 
mdq m dp^ 

T BK --knT—— (IRc) 

L * " kBl dp % (18C) 

LT = -ksT^^- (18d) 
dq op^ 

jBK = J_dGidKid_ _ 

5 dp dpc_ d( dp^ dp 

G 2 dK 2 d i dG 2 dK 2 d 



BK 



dp^ dq dq dp^ d£ 

Upon choosing a symmetric Trotter factorization of the BL Liouville operator based on the 
decomposition 

£ BK = X>* K (19) 



a=l 



a measure-preserving algorithm can be produced in full generality. 

In practice, a choice of K X) K 2 , G\, G 2 must be made in order obtain explicit formulas. 
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In this paper, we make the following simple choices: 



G l =p 
G 2 = q 

p\ 



K 2 



P\ 
2 



In terms of Eqs fl20aH20dl) . the antisymmetric BK tensor becomes 



B BK 



and the Hamiltonian reads 

H BK 



1 -q 
1 
1 
-10 -p 
-1 p 
q 0-100 



2 2 

H(q,p) + -p- + p-+k B T(C + Z) ■ 



2m^ 2m^ 



The split Liouville operators now simplify as follows: 



' BK 



: BK 



: BK 



: BK 



L 



BK 



_dV d dV d 
dq dp ^ dq dp^ 
p d p 2 d 
mdq m dp^ 

-k B T— 

-k B T— 

% 

Pc, d Pc, d 
m^d( dp 
Bt 9 vc d 



dp v 
p\ d 



(20a) 
(20b) 

(20c) 
(20d) 



(21) 



(22) 



(23a) 
(23b) 
(23c) 
(23d) 

(23e) 
(23f) 



dq <9£ dp x 

For the purposes of defining an efficient integration algorithm, we combine commuting 
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Liouville operators as follows: 



where 



Defining 



j BK 



-BK 



j BK 

C = 



Lf K + L? 
dp 



d 



BK 

3 



i- |- F 

m dq Pc dpc 

f BK , r BK 
L h + ^6 

Pc d pe d pt d pc d 
^TT + 7^7 + 



F(q) = -dV/dq 



dq 



PC 



U* K (r) 



P 7 rp 

k B T . 

m 



exp 



tL 



BK 



(24a) 



(24b) 



(24c) 



(25a) 
(25b) 

(25c) 
(26) 



where a = A,B,C, one possible reversible measure-preserving integration algorithm for the 
BK thermostat is then 



U(r 



,BK 



' B 



x UT (t) 



iK bk (iK© 



- d K o c c bk © d K v , 



(27) 



Using the so-called direct translation technique |36j we can expand the above symmetric 
break-up of the Liouville operator into a pseudo-code form, ready to be implemented on the 
computer: 



q -> q + $ 
PC -> Pc + \ F pc 



uE K (I) 
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p — > pexp 
q — > gexp 



T Pj 

2 

T Pg 

2 



^ ^ T 2 m c 
£ £ + r PL 



ttBK 



(I) 



9 -> g+ IiL 



4 TO 



->• R + i f p 



= ^ K (i) 



p — >■ p + tF 
P£ ->■ P£ + 



^ ^ 4 m 

' I 4 rrirj 



Pc ->• pc + 



PC 



(5) 



p — > pexp 
g — >■ gexp 

C -x C + | 



£ £ + IPL 



r PC 
2 

_I PL 

2 



PC 



PL 



^ (1) 



i ^ 4 m 

pc ->■ pc + i f p 



■ /r BK 



IV. BULGAC-KUSNEZOV-NOSE-HOOVER DYNAMICS 



The BKNH Hamiltonian 



tf BKNH = H{q,p) + + + A. + kB T(C + + 2^ 



(28) 



is simply the BK Hamiltonian augmented by the Nose variables {r),p v ) with mass m v . With 
the antisymmetric BKNH tensor 



B 



BKNH 















1 





-G 2 




















dGi 
dp 


























8G 2 
8q 


























1 


-1 














-Gt 











dGi 
dp 








G a 








-p< 


G 2 





dG 2 

dq 














-Pt 











-1 





PC 


PC 






(29) 



we obtain from Eq. ([3]) equations of motion for the phase space variables x 
(q,(,tV,P,P0Pt,Pv) = (xi,x 2 ,x 3 ,X4,x 5 ,x 6 ,x 7 ,x 8 ): 

OH G 2 (q,p)dK 2 



C 

i 
q 

p 

PC 

p\ 

Pv 



dp dp^ 

1 8G X 8K X 
rriQ dp dp^ 

1 dG 2 dK 2 
dq dp^ 

Pv 
m v 

dH G 1 (q,p)dK 1 



G l 



G 



dq 
dH 
dp 
dH 



dp^ 

dGi p v 

PC, 



dq 
P( dKi 



k B T 
k B T 



dp 

dG 2 p v 

pc 



dq 
P( dK 2 



2k B T 



(30a) 
(30b) 
(30c) 
(30d) 

(30e) 
(30f) 

(30g) 
(30h) 



dp^ dp^ 

Here, a single Nose variable is coupled to both of the BK demons ( and £. The associated 
invariant measure is discussed in Appendix O 
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A. Algorithm for BKNH dynamics 



The Hamiltonian can be split as 

Hf KNH = V(q) (31a) 

^BKNH = f_ (31b) 

Hf KNH = k B T( (31c) 

Hf KNK = k B Ti (31d) 

^bknh = Kiipd (31e) 

H BKNR = KltPd (31f) 
2 

rrBKNH _ Py /oi _\ 

H " -2m~ v (31g) 

# 8 BKNH = 2k B Tr] (31h) 



The measure-preserving splitting 19j of the Liouville operator 

a rxBKNH a 

,on n a 



yields 



l bknh = + Q W_ A ( 33a ) 

1 <9g dp dq dp^ 

L BKKK = ^d +G ^d 

mdq m dp^ 

L ™ = (33c) 

l bknh = ^_ _ ChdKid_ + p^dKi_d_ 

dp dp^ <9£ ttiq dpc_ dp dpc_ dp n 

g 2 8k 2 d i aG 2 ax 2 a P€ dK 2 d 

H ~ 7\ — 7T7 H — q — 7^— (331) 



BKNH 



mg <9p^ <9g dq dp^ <9£ dp^ dp n 



L bknh = P^d _ p^ d _ p^ d 
m v orj dp^ m v dp^ 

L™ = -2k B T^- . (33h) 

OPri 
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At this stage, we leave the general formulation and adopt the particular choice of K\, 
K 2 , Gi, and G 2 given in Eq. (TSUI) . The antisymmetric BKNH tensor becomes 



^BKNH 















1 





-q 




















1 


























1 


























1 


-1 














-p 











-1 








p 








-p< 


<1 





-1 














-Pi 











-1 





p< 








(34) 



and the Hamiltonian simplifies to 



H 



BKNH 



H{q,p) + 



pI 



+ 



2m^ 2m, 



P 2 V 
2m r 



+ k B T(( + + 2 B Tr/ . 



(35) 



The split Liouville operators are now 



LI 



BKNH 



L 



BKNH 



_dV d_ dV d 
dq dp ^ dq dp^ 
p d ^p 2 d 



BKNH 



; BKNH 



: BKNH 



: BKNH 



^BKNH f'l " _ ni " ft 



mdq m dp^ 

-k B T— 

dP( 

-k B T— 

opt 

_ el ll H 

m^ d( dp 

Pi d p£ d 
-Q— + 



P^_d_ 

rriQ dp v 
P\ d 



+ 



mg dq m^ d£ m^ dp v 

p v d p v d p v 
PC 



L 



BKNH 



rriri dr] 
-2k B T 



m 
d 

dp v 



d 



dp ( 



m. 



-Pi 



(36a) 
(36b) 
(36c) 
(36d) 
(36e) 
(36f) 
(36g) 
(36h) 



For the purposes of defining an efficient integration algorithm, we combine commuting Li- 
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ouville operators as follows: 



^BKNH _ ^BKNH _|_ j^BKNH _|_ ^BKNH 



dp m v or] m, n dp^ 



j BKNH _ j BKNH , f BKNH 



BKNH 



J 2 

p d 







J c 



mdq^ Pc ~dp^ 

j BKNH | j BKNH , f BKNH 

_hl — L\ + El A + elIL _|_ p — 

dp dq d( d£ Pv p v 



where 



dV 
dq 
dV 



Pn 



P j rp 

k B T 

m 



In La there appears an operator with the form 

Pk 



m k * J dp. 



d 



(37a) 
(37b) 
(37c) 

(38a) 
(38b) 
(38c) 
(38d) 

(39) 



where (k,i) = (%, £) f° r La- The action of the propagator associated with this operator on 
Pi is derived in Appendix |XJ and is given by 

(40) 



T - T £k- T Pk 

e 1 'Pi = Pie m k + rF Pi e 2m '- 



k t 



Pk 

2mi 



sinh 



, Pk 
2mi 



The apparently singular function 



. Pk 
2mi 



sinh 



. Pk 
2rrii 



(41) 



is in fact well behaved as p k — > 0, and can be expanded in a Maclaurin series to suitably 
high order {3?]]. In our implementation we used an eighth order expansion. 
The propagators for the BKNH dynamics can now be defined as 



U* KNK (t) = exp 



tL 



BKNH 



(42) 
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where a = A, B, C. One possible reversible measure-preserving integration algorithm for 
the BKNH thermostat can then be derived from the following Trotter factorization: 



C/(r)™ = f/™ (I) UE K ™ (I) f/I KNH (£) 
x C/™(r) 

_ rrBKNH / TN \ rrBKNH / T ^ rrBKNH / T ^ 

x U B {-) U c {-) U B {-) 
The direct translation technique gives the following pseudo-code: 

(5) 



9 -> 9+ i-^ 



pc ^ pc + i f pc 



rrBKNH /> N 



(43) 



p — )■ pexp 
q — >• gexp 



T PC 
2 TO,J 

T Pg 

2 m,£ 



Pr, -> Pt, + i-PJ 



P< 



rrBKNH 



(5) 



9 -> 9+ 5* 
PC ^ PC + \ F i 



PC 



jyBKNH (t 



(i) 



p — r p + rF(q) 
Pa ->■ PC + 



PC — ^ Pc ex P 



. P>7 



> : C/™( 



9 

Pc ^ Pc + i F P C 



rrBKNH (t 



(5) 



p — > pexp 
q — > gexp 



T PC 
2 T7l£ 

2 



C -r <f + 

^ ^ <» T 2 m f 
Pr; ^ Pr; ~r" 2 F Pv 



(i) 
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H 1 4 m ^ . rr BKNH 



^BKNH (i) 



V. BULGAC-KUSNEZOV-NOSE-HOOVER CHAIN 



For simplicity, we explicitly treat only the case in which the pq an d p$ demons are each 
coupled to a standard NH thermostat (length one). It would be straightforward to couple 



each of the demons to NH chains 
follows. Define the Hamiltonian 



32] , and the general case can be easily inferred from what 



H 



BKNHC 



Px 



+ -f2^ + k B T(( + Z + r ] + x) ■ (44) 



2m„ 2m 



m^ m^ 

Upon defining the phase space point x = (q, £, £, rj, XiPiP^Pi-iPniPx) 
(xi, x 2 , x 3 , x 4 , x 5 , x e , x 7 , Xs, x 9 , Xiq) and the antisymmetric BKNHC tensor 



B 



BKNHC 


















1 





-G 2 


























dGi 
dp 
































dG 2 
dq 
































1 
































1 


-1 

















-G, 














_8Gi 

dp 











G 1 








-PC 





G 2 





8G 2 
dq 




















-Pi 











-1 








PC 























-1 








Pi 









associated non-Hamiltonian equations of motion are 

^^BKNHC 



X i 



BKNHC 



dx; 



(45) 



(46) 



with i 



10. 
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A. Algorithm for BKNHC chain dynamics 



Splitting the BKNHC chain Hamiltonian as 

^bknhc _ 



^jBKNHC 



rrBKNHC 
h 3 



H 



BKNHC 



BKNHC 



H 



rrBKNHC 
H 6 



H- 



BKNHC 



H. 



BKNHC 



V(q) 

2m 
k B T( 

k B Tti 
KM) 

P\ 

2m v 
k B Trj 



H, 



BKNHC 



P 



x 



2m, 



H 



BKNHC 
10 



"X 



we obtain the corresponding measure-preserving splitting of the Liouville operator 

^BKNHC q 



10BKNHC 1 

ij 



dXn 



dxi 



(47a) 

(47b) 

(47c) 
(47d) 

(47e) 
(47f) 

(47g) 
(47h) 
(47i) 
(47j) 

(48) 



At this stage we go directly to Eqs ( 12U|) . The antisymmetric Nose-Hoover-Bulgac- 
Kusnezov tensor becomes 



^BKNHC 


















1 





-q 


























1 
































1 
































1 
































1 


-1 

















-p 














-1 











p 








-Pc 





Q 





-1 




















-Pi 











-1 








PC 























-1 








Pi 









(49) 
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the Hamiltonian 

tf™c = H (q ,p) + A. + A. + ^ + A- + k B T( C + f +„ + 

and associated Liouville operators 

r bknhc _ dV d dV d 

L l ~ o "o + 9" 



^BKNHC _ J- " | 



<9g dp <9g <9p^ 



m dq m dp^ 
3 — —Kb-! 



^BKNHC = u T 



dp c 

4 — —kb-< 



L BKNHC = T 



p^ d Pc, d p\ d 



L BKNHC = ^_ _ ^ ~ + 



P(_ d p^ d p\ d 



^BKNHC _ _£^_gj_ _|_ _|_ 

6 rri£ dq d£ dp x 

^bknhc _ Pv_j^_ HlL v - 

7 m v Or] m v dp^ 

^BKNHC = u T 



dp v 



^BKNHC _ P*_Jl Px 9 

9 m x dx rn x dp^ 



Z BKNHC = _ kBT . 



d 



We combine commuting Liouville operators as follows: 



r BKNHC _ f BKNHC 



t BKNHC 



_|_ j^BKNHC _|_ j^BKNHC 




r BKNHC _ 
L c = 



. ^BKNHC _|_ j^BKNHC _|_ j^BKNHC 
p_8_ _P^^_ /_ Prj_ 

m dq m v di] \ m v 




j BKNHC | j BKNHC , f BKNHC 

p c Pz d p c d 

Pa IJT + ^7 ' 

op oq oQ 



f BKNHC 
" ^10 

Pi d 



m £, °q p v 



d 
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where 



F(q) = - 



dV_ 

dq 
dV 



Pi 



pq 



Pn 



Px 



q^. — k B T 
oq 



knT 



m 



knT . 



(53a) 
(53b) 
(53c) 
(53d) 

(53e) 



Both in L^ KN and L^ KNHC there appears an operator with the form 



(54) 



where (k,i) = (x, £) for La and (k,i) = (r],() fo r Lb- Again following the derivation in 
Appendix [Aj we find 

-i 



e rLl pi = pie ' m k + rF Pt e ' 2m fc I r 



Pk 
2rrii 



sinh 



, Pk 
2mi 



(55) 



The function ( ) s ^ nn 



. Pk 

2m k 



is treated through an eighth order expansion 



37] . 



The propagators 



U, 



BKNHC . 



T) = exp 



tL 



BKNHC 



(56) 



with a = A,B,C can now be introduced. One possible reversible measure-preserving inte- 
gration algorithm for the BKNHC chain thermostat is then 



rj/ r \BKNHC _ ryBKNHC f ~[\ ryBKNHC 

\ ) b y A j c \ 2 

x U* KNHC (t) 



v r /-BKNHC f T \ rrBKNHC f T \ rrBKNHC f T ^ 

X UJ Uc \2J U, 



(57) 



In pseudo-code form, we have the resulting integration algorithm: 



q -> q + j- 
v ^ v + j— 

i i 1 4 m- 



PC -> Pc e 



1 

t Py t Py 

4 m„ _|_ L J? o 4 2m„ I T Py 



V, . ( rBKNHl r 



+ lF Pi e l2 "M^) sinh 



4 2m r 



T Py 
4 2m, 



(9 
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p — > pexp 
q — > q exp 



T PC 

2 
2 



Px ^ Px + \ F P X 



> : f/^KNHC (z) 



1 n 1 4 m 



77 — >• 77 + 



T P?7 

4 m,, 

r Pt? t Py 

4 m, I T J? p 4 2m„ 

4 X >C 



p — > p + rF 
y — >• y + r-^- 

^ ^ m x 



4 2m r , 



sinh 



T Py 
4 2m„ 



sinh 



T 



Px 
2m» 



T PT? / \ 

4 2m, [ T P»7 \ 

\4 2m,,/ 



-1 



sinh 



T Py 
4 2m, 



p — > pexp 
g — >■ gexp 

P v ~> P v + \F t 
Px ^ Px + l F i 



T PC 

2 m,Q 

T P( 

2 m^ 



Py 
Px 



9 


->■ 


9 + 


11 

4 m 


• 77 


->■ 


7/ + 


l_PiL 

4 m, 

T PT) 


p< 


->■ 




4 m, 



4^PC ( 



T Pi? / \ — 1 

4 2m, [ r Py \ 
\i 2m v J 



L 

sinh 



T_Py_ 

4 2m, 
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VI. NUMERICAL RESULTS 



In its simplicity, the dynamics of a harmonic mode in one dimension is a paradigmatic 
example for checking the chaotic (ergodic) properties of constant-temperature phase space 
flows and the correct sampling of the canonical distribution. It is well known that it is 
necessary to genera . 
Nose-Hoover chain 
systems such as the harmonic oscillator. 

Some time ago, BK dynamics was devised to provide a deterministic thermostat for 



ize basic Nose- Hoover dynamics Is), 9] to thermostats such as the 
1, [3?! in order to produce correct constant-temperature averages for 



systems such as classical spins 



23 



241 ] . To ensure efficient thermostatting, BK found it 



necessary to introduce several 'demons' per thermostatted degree of freedom, where each 
demon was taken to have a different and in principle complicated coupling to the system 



degree of freedom 



23l . |24J. In the present work, we keep the form of the system-thermostat 



coupling as simple as possible, in order to facilitate the formulation of explicit, reversible 
and measure-preserving integrators [3]. It is then of interest to investigate the ability of our 
BK-type thermostats to produce the correct canonical sampling in the case of the harmonic 
oscillator. Interest in harmonic modes is also justified by the possibility of devising models 
of condensed matter systems in terms of coupled spins and harmonic modes, as already done 
in quantum dynamics with so-called spin-boson models [38?]. We therefore investigate the 
performance of our integration schemes on the simple one-dimensional harmonic oscillator. 

For the particular calculations reported here, the oscillator angular frequency, all masses 
and ksT were taken to be unity. The time step in all cases was r = 0.0025, and all runs were 
calculated for 10 6 time steps, starting from the same initial conditions: harmonic oscillator 
coordinate q = 0.3, all other phase space variables zero at t — 0. 

The measure-preserving algorithms derived here result in stable numerical integration for 
all the three cases treated: BK, BKNH, and BKNHC chain dynamics. Figure [T] shows the 
three extended Hamiltonians (normalized by their respective initial time value) versus time. 
All three Hamiltonians are numerically conserved by the corresponding measure-preserving 
algorithm to very high accuracy (which is maintained in all the three cases). 

However, the basic BK phase space flow is not capable of producing the correct canonical 
sampling for a harmonic mode. This can be easily checked since the canonical distribution 
function of the harmonic oscillator is isotropic in phase space and its radial dependence 



23 



can be calculated exactly. Details of this way of visualizing the phase space sampling have 



already been given in 



14 



151 ] . Figure [U displaying the comparison between the theoretical 



and the calculated value of the radial probability in phase space, clearly shows that the 
BK dynamics is not able to produce canonical sampling. A look at the inset of Fig. [2, 
showing the phase space distribution for the harmonic mode, also immediately shows that 
the dynamics is not ergodic. 

The same analysis has been carried out for BKNH and BKNHC phase space flows, and 
these are displayed in Fig [3] and Fig HJ respectively. Within numerical errors, both BKNH 
and BKNHC thermostats are able to produce the correct canonical distribution for the stiff 
harmonic modes. 

Introduction of a single, global Nose-type variable in the BKNH thermostat effectively 
introduces additional coupling between the two demon variables. The effectiveness of the 
BKNH thermostat is consistent with our findings (results not discussed here) that introduc- 
tion of explicit coupling between demons in BK thermostat dynamics also leads to efficient 
thermostatting of the harmonic oscillator. 
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VII. CONCLUSIONS 



We have formulated Bulgac-Kusnezov 23|, |24|, Nose-Hoover controlled Bulgac-Kusnezov, 



and Bulgac-Kusnezov-Nose-Hoover chain thermostats in phase space by means of non- 



Hamiltonian brackets 
rithms 



14 



151 ] . We have derived time- reversible measure-preserving algo- 



191 ] for these three cases and showed that additional control by a single Nose-Hoover 



thermostat or independent Nose-Hoover thermostats is necessary to produce the correct 
canonical distribution for a stiff harmonic mode. 

Measure-preserving dynamics of the kind discussed here is associated with equilibrium 
simulations (where, for example, there is a single temperature parameter T). Stationary 
phase space distributions associated with non-equilibrium situations are much more com- 



plicated than the smooth equilbrium densities analyzed in the present paper [111, [39|, 140 1. 
Nonequilibrium simulations of heat flow could be carried out by extending the present ap- 
proach to multimode systems (e.g., a chain of oscillators) coupled to BK-type demons with 



41 



43]. 



associated NH thermostats corresponding to two different temperatures 

The techniques presented here for derivation and implementation of thermostats have 
been shown to be efficient and versatile. We anticipate that analogous approaches can be 
usefully applied to systems of classical spins coupled to both harmonic and anharmonic 
modes. 
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Appendix A: Operator formula 



We wish to determine the action of the propagator associated with the Liouville operator 
Eq. (|39|) . This is equivalent to solving the evolution equation (recall i ^ k) 



dpi ( Pk . „ 

— — = Pi + iy 

at V m k 



from t = to t=r. Integrating, we have 



m k p k 

In pi + F Pi 

Pk \ m k 



giving 



Pi{T) = exp 



Pk 



Pi 

-rp k /m k ^ 



sinh 



P 



e -rp k /m k + r _p T 2^ 



T 



2mt 



T 



Pfc 
2m k 



Appendix B: Invariant Measure of the BK phase space flows 



(Al) 

(A2) 

(A3a) 
(A3b) 

(A3c) 



The phase space compressibility of the phase space BK thermostat is 



dBffdH BK 1 dddKi 1 dG 2 dK 2 



dxi dxi 
Upon introducing the function 



dp dpQ dq dp^ 



H, 



BK 



H 



Ki + K2 



one can easily find that 



^BK 



1 dH* K 



ksT dt 

so that the invariant measure in phase space reads 



dfi = dx exp 



dtKftK 



= dx exp [-/3#| K ] 

= dx exp[-/3H BK ] exp[C + f] 



(Bl) 

(B2) 
(B3) 

(B4a) 

(B4b) 
(B4c) 



as desired. 
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Appendix C: Invariant Measure of the BKNH phase space flows 



The phase space compressibility of the NH controlled Bulgac-Kusnezov thermostat is 

<9ffg KNH dH BKNH = _J_dG 1 dK 1 _ ^OGidKz _ ^ 
dxi dxi dp dp^ dq dp^ m v 

Upon introducing the function 



if^ KNH = H + — + — + (C2) 

we have 

KBKNU = J^f— (C3) 



so that the invariant measure in phase space is 



dfi = dx exp 



— / dtK 



BKNH 



we have 



so that the invariant measure in phase space reads 



dp, = dx exp 



— / dtKBKNHC 



(C4a) 



= dx exp [-f3H* KNR ] (C4b) 
= dx exp[-/3# BKNH ] exp[C + £ + 2??]. (C4c) 

Appendix D: Invariant Measure of the BKNHC chain phase space flows 

The phase space compressibility of the Nose-Hoover-Bulgac-Kusnezov chain is 
9£§ knhc c>#bknhc 1 dddK, 1 OG 2 dK 2 Pv Px 

^BKNHC — ^ ^ — ^ ^ K ~ {L>1) 

oxi oxi op op^ oq op^ m v m x 

Upon introducing the function 



= H + — + — + + (D2) 
2m v 2m x 



1 dH* K 

^BKNHC = j-f— (D3) 



(D4a) 



dx exp [-/3#™ c ] (D4b) 
da; exp[-/3# BKNHC ] exp[C + £ + V + x\- (D4c) 
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FIG. 1: Comparison of the total extended Hamiltonian versus time (normalized with respect to its 
value at t = 0) for the harmonic oscillator undergoing simple Bulgac-Kusnezov dynamics (H BK ), 
NH controlled Bulgac-Kusnezov dynamics (i7 BKNH ), and Bulgac-Kusnezov-Nose-Hoover chain dy- 
namics curves have been displaced vertically for clarity. The time-reversible 
measure-preserving algorithms developed in this paper conserve the extended Hamiltonian to high 
accuracy in all three cases. 



30 




0.5 1 1.5 2 2.5 3 3.5 4 

r 

FIG. 2: Radial phase space probability for a harmonic oscillator under Bulgac-Kusnezov dynamics. 
The continuous line shows the theoretical value while the black bullets display the numerical results. 
The inset displays a plot of the phase space distribution of points along the single trajectory used 
to compute the radial probability. 
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FIG. 3: Radial phase space probability for a harmonic oscillator under Nose-Hoover controlled 
Bulgac-Kusnezov dynamics. The continuous line shows the theoretical value while the black bullets 
display the numerical results. The inset displays a plot of the phase space distribution of points 
along the single trajectory used to compute the radial probability. 
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FIG. 4: Radial phase space probability for a harmonic oscillator under Bulgac-Kusnezov-Nose- 
Hoover chain dynamics. The continuous line shows the theoretical value while the black bullets 
display the numerical results. The inset displays a plot of the phase space distribution of points 
along the single trajectory used to compute the radial probability. 
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